clear; close all;
saveArg = 1;


%% Plotting

medSize = 70;
w = 3;
quivSkip = 7;
quivSize = 5;

labelSzAx = 16;
labelSzXSec = 18;

MadMedTheshold = 0.6;

% Circulation
data = importdata('Exp_21OCTSolid\Data_03-15000-15100\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);

thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);

indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

figure(1000); clf;
set(gcf,'units','normalized','OuterPosition',[0.05 0.1 0.9 0.6]);
subplot(1,3,1);
scatter(xMed,yMed,30,thetaMed,'filled'); hold on;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),quivSize,'k');
set(gca,'Fontsize',labelSzAx);
axis equal;
xlim([-100,150]);
ylim([-300,-50]);
set(gca,'XTick',-100:100:100);
set(gca,'YTick',-300:100:-50);
ylabel('y (mm)');
caxis([0,pi]);
tx = text(-80,-70,'a'); set(tx,'Fontsize',26);



quivSkip = 9;
quivSize = 3;

% Buoyant
data = importdata('Exp_Prev_C1\Data_24-10200-10300\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);



thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);

indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
end
for i = 1:length(indY)
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

subplot(1,3,2);
scatter(xMed,yMed,30,thetaMed,'filled'); hold on;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),5,'k');
set(gca,'Fontsize',labelSzAx);
xlabel('x (mm)');
axis equal;
xlim([-200,200]);
ylim([-300,100]);
set(gca,'XTick',-200:200:200);
set(gca,'YTick',-300:200:100);
caxis([0,pi]);
tx = text(-165,68,'b'); set(tx,'Fontsize',26);



% Feeder
data = importdata('Exp_Prev_C1\Data_27-11400-11500\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);



thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);


indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
end
for i = 1:length(indY)
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

subplot(1,3,3);
scatter(xMed,yMed,30,thetaMed,'filled'); hold on;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),5,'k');
set(gca,'Fontsize',labelSzAx);
axis equal;
ylim([-300,100]);
xlim([-200,200]);
set(gca,'XTick',-200:200:200);
set(gca,'YTick',-300:200:100);
caxis([0,pi]);
tx = text(-170,70,'c'); set(tx,'Fontsize',26);
c=colorbar;
ylabel(c,'\theta (deg)');
set(c,'Ticks',[0,pi],'TickLabels',[0,180]);

pause(0.1);
x1 = 0.08; x2 = 0.37; x3 = 0.66;
y1 = 0.18;
wd = 0.21; ht = 0.73;
subplot(1,3,1); set(gca,'Position',[x1 y1 wd ht]);
subplot(1,3,2); set(gca,'Position',[x2 y1 wd ht]);
subplot(1,3,3); set(gca,'Position',[x3 y1 wd ht]);
set(c,'Position',[0.91 y1 0.02 ht]);




%% Saving

if saveArg == 1
    outFile = 'results_plot';
    print(gcf,outFile,'-djpeg','-r300');
    saveas(gcf,outFile,'fig');
    saveas(gcf,outFile,'epsc');
end

